##Replication File for "Authoritarian Audiences, Rhetoric, and Propaganda in International Crises: Evidence from China"
##Authors: Jessica Chen Weiss and Allan Dafoe

##Replication Hypothetical History


#Packages
library(ggplot2)
require(reshape2)
require(coefplot)
library(robustbase)

#hyp.Rdata/real.Rdata cannot be shared due to privacy concerns, instead the de-identified datasets with simulated demographics (hyp_toshare and real_toshare) can be used to replicate the paper's results.

setwd("~/Dropbox/Dafoe-Weiss-Empirics/18-04-12-ISQ-Raluca/Code/Replication File ISQ/")

#------------------------------------------------------------------------------------------
## Constructing Publishable dataset. 
#------------------------------------------------------------------------------------------
#load("./MASTER CODE/hyp.Rdata")
#keepvars <- c("asc", "asc0", "pro", "prot", "com", "mob", "eli.f", "eli.c", 
#    "authoritarian", "ally", "capabilities", "salience", "his",
#    "pre.questions", "asc.or", "asc0.v2", "na1.v2", "na2.v2", "na3.v2",
#    "na2.v.dn", "na3.v.dn", "start.time", "start.time.n", "start.time.n2", "start.time.n3",
#    "start.time.swd", "att", "rac", "rac.or", "ra3c", "na1.v", "na2.v", "na3.v", "a.pcono") ##variables that are not identifiable


#data.share <- d[, names(d) %in% keepvars]
#write.csv(data.share, "./MASTER CODE/hyp_toshare.csv") 


d <- read.csv("./hyp_toshare.csv") #analysis data keeps only non-identifiable variables. 

#generating random numbers of identifiable variables
set.seed(2018)

d$gender <- sample(c(0,1), size=nrow(d), replace=TRUE, prob=c(.39,.61))
d$educ <- sample(seq(from = 1, to = 7, by = 1), size = nrow(d), replace = TRUE, prob=c(0.002, 0.00, 0.402,  0.067, 0.470, 0.054, 0.005))
d$age <- sample(c( 6,11,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32, 
                   33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50, 
                   51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68 ,
                   69,70,71,72,73,85),
                size=nrow(d), replace=TRUE, 
                prob=c(0.0003, 0.0003, 0.0006, 0.0025, 0.0039, 0.0066, 0.0063, 0.0132, 0.0113, 0.0121, 0.0215, 0.0157, 0.0212, 
                       0.0251, 0.0240, 0.4307, 0.0196, 0.0317, 0.0287, 0.0259, 0.0389, 0.0245, 0.0187, 0.0207, 0.0204, 0.0185,
                       0.0143, 0.0149, 0.0110, 0.0105, 0.0141, 0.0094, 0.0066, 0.0063, 0.0052, 0.0050, 0.0036, 0.0091, 0.0066,
                       0.0044, 0.0033, 0.0025, 0.0033, 0.0025, 0.0039, 0.0025, 0.0028, 0.0022, 0.0017, 0.0017, 0.0022, 0.0011,
                       0.0014, 0.0011, 0.0003, 0.0014, 0.0014, 0.0003, 0.0006, 0.0003))
d$partner <- sample(c("A", "B"), size=nrow(d), replace=TRUE, prob=c(.551, 0.449)) 


#Gen other variables 
d$gender.o <- d$gender
d$gender.m <- is.na(d$gender)
d$gender[d$gender.m==TRUE] <- 1

d$age.o <- d$age
d$age.m <- is.na(d$age)
d$age[d$age.m==TRUE] <- 30

d$educ.o <- d$educ
d$educ.m <- is.na(d$educ)
d$educ[d$educ.m==TRUE] <- 3


write.csv(d, "./MASTER CODE/hypanalysis.csv")
save(d, file="./MASTER CODE/hypanalysis.Rdata")
  

#-----------------------------------------------------------------------------------------
#Replication Main Text and Appendix
#-----------------------------------------------------------------------------------------


{{format(Sys.Date(), format="%B %d %Y")}}

#+ opts, include=FALSE,eval=TRUE
knitr::opts_chunk$set(eval=TRUE, echo=TRUE, error=TRUE,message=FALSE, warning=FALSE,
                      fig.width=15, fig.height=5,width=60,dev="pdf"
)


rm(list = ls())

setwd("~/Dropbox/Dafoe-Weiss-Empirics/18-04-12-ISQ-Raluca/Code/Replication File ISQ/")

load("./MASTER CODE/hypanalysis.Rdata")
rb <- rep(1, length(d[,1]))



#Figures location
fig.loc <- "../../Figures/"

#Packages
require(ggplot2)
require(reshape2)
require(coefplot)
require(robustbase)
require(rms)
require(sandwich)
require(stargazer)
w <- 5
h <- 4

#finding value of Z that gives 10% one sided test
xs <- seq(-2,0, 0.01)
xs[which(pnorm(xs)-0.05==min(abs(pnorm(xs)-0.05)))]
pnorm(-1.64)
pnorm(-1.96)

innerCI.n <- 1.64 #set so exclusion of CI implies two sided significance of 0.1
outerCI.n <- 1.96 #set so exclusion of CI implies two sided significance of 0.05




# Customized Stargazer Function, without asterisk inflation
stargazerAD <- function(table, title="Default", filename="Figures/default.tex", dep.var.labels=NULL) {
  a <- stargazer(table, title=title,  out=filename, 
                 align=TRUE, star.cutoffs = c(0.1, 0.05, 0.01, 0.001), star.char = c("\\dagger", "*", "**", "***"), 
                 notes="$^{\\dagger} p<0.1$; $^{*} p<0.05$; $^{**} p<0.01$; $^{***} p<0.001$", notes.append=FALSE, single.row=TRUE)
  return(a)
}



results.fun <- function(data, r, model, title.w, fig.loc, coefficients=NULL, innerCI.n=1.64, outerCI.n=1.96, w=10, h=6) {
  d <- data
  #rob.results <- robcov(ols(model, dat=d[r,], x=TRUE))
  m1 <-  lm(model,  dat=d[r,])
  homo.results <- summary(m1)
  n <- length(d[r,1])
  title <- paste(title.w, ", n=", n, sep="")
  pf <- coefplot(m1, title=title, 
                 intercept=FALSE, lwdOuter=0.5, innerCI=innerCI.n, outerCI= outerCI.n,
                 coefficients=coefficients) + theme_bw()
  fig.results <- pf
  list.results <- list("homo.results"=homo.results, "title"=title, "fig.results"=fig.results, "model"=m1, m1)
  return(list.results)
}



## Number of observations ##
n.obs <- sum(table(d$asc))
cat(n.obs, file=paste(fig.loc,"Hyp-n.tex", sep=""))

n.obs.att <- sum(d$att) 
cat(n.obs.att, file=paste(fig.loc,"Hyp-n-att.tex", sep=""))



#From PAP: 
#For our primary specification we will estimate the effects of our manipulations using OLS with heteroscedasticity consis- tent standard errors, conditioning on the other experimental conditions to increase power. Hypothesis tests will be one-sided when appropriate. We will also report our estimates of the average causal effect when not conditioning on anything else, an estimand termed the average marginal component effect by Hainmueller, Hopkins, and Yamamoto.37
#Our primary test of these hypotheses will involve OLS regression, controlling for the other experimental conditions, and in the real-history design for the pre-scenario answers to the same questions. 
#robcov(ols(y ~ a.1 * a.2, toyData, x = TRUE)) ## heteroskedastic (ATE)

#' ### Results for Approval ###

#r <-  rep(1, length(d[,1]))
r <- !is.na(d$asc) & rb

type <- "Hyp"
outcome <- "Approval"
cont <- ""
title.w <- paste(type, ", Effect on ", outcome, cont, sep="")
title.f <- paste(type, "-", outcome, cont, sep="")

coefficients.list<- c("his", "pro", "com", "mob", "eli.f", "eli.c")
model <- asc ~ pro + prot + com + mob + eli.f + eli.c + 
  authoritarian + ally + capabilities + salience + his +
  pre.questions + asc.or 
results <- results.fun(data=d, r=r, model=model, title.w=title.w, fig.loc=fig.loc, coefficients=coefficients.list, innerCI.n=innerCI.n, outerCI.n=outerCI.n, w=w, h=h)
results[1]
results[2]
results$fig.results
coefplot(lm(model, data=d[r,]), title=title.w, 
         intercept=FALSE, lwdOuter=0.5, innerCI=innerCI.n, outerCI= outerCI.n, data=d) + theme_bw()

pdf(paste(fig.loc,title.f,".pdf", sep=""), width=w, height=h)
results$fig.results
dev.off()
m11 <- results[4]

pro.results <- coef(summary(results$model))["pro",]
toprint <- paste("p_{+}=",round(pro.results[4]/2,3),  sep="")
write(toprint, file=paste(fig.loc,type,outcome,"p",".tex", sep=""))

statement.results <- coef(summary(results$model))[c("com", "mob", "his", "eli.f","eli.c"),]


#Checking if provocation or statement of commitment interacts with eli.f or eli.c.
#Audience costs arise if pro or com decrease approval (negative effect). Eli would reduce audience costs if they have a positive interaction. 

model <- asc ~ com + pro + prot + com*eli.f + com*eli.c + pro*eli.f + pro*eli.c + mob + eli.f + eli.c +
  authoritarian + ally + capabilities + salience + his +
  pre.questions + asc.or 
coefplot(lm(model, data=d[r,]), title=title.w, 
         intercept=FALSE, lwdOuter=0.5, innerCI=innerCI.n, outerCI= outerCI.n, data=d) + theme_bw()

#Result: largely insignificant, but if anything statement of commitments and elite cue framing negatively interact. So people disapprove more if you make a statement of commitment, and then offer a biding time argument for justifying backing down, perhaps because you are revealing inconsistency (but why not same result for eli.c?). 



r <- !is.na(d$asc) & rb

type <- "Hyp"
outcome <- "Approval"
cont <- "Covariates"
title.w <- paste(type, ", Effect on ", outcome, ", Controlling ", cont, sep="")
title.f <- paste(type, "-", outcome, cont, sep="")

#coefficients.list<- c("his", "pro", "prot", "com", "mob", "eli.f", "eli.c")
coefficients.list<- c("his", "pro", "com", "mob", "eli.f", "eli.c")
model <- asc ~ pro + prot + com + mob + eli.f + eli.c + 
  authoritarian + ally + capabilities + salience + his +
  pre.questions + asc.or +  partner + asc.or + asc0.v2  +
  na1.v2 + na2.v2 + na3.v2 + na2.v.dn + na3.v.dn + gender + educ + age + age.m + gender.m + educ.m + 
  start.time.n + start.time.n2 + start.time.n3 + start.time.swd 
results <- results.fun(data=d, r=r, model=model, title.w=title.w, fig.loc=fig.loc, coefficients=coefficients.list, innerCI.n=innerCI.n, outerCI.n=outerCI.n, w=w, h=h)
results[1]
results[2] 
results$fig.results
coefplot(lm(model, data=d[r,]), title=title.w, 
         intercept=FALSE, lwdOuter=0.5, innerCI=innerCI.n, outerCI= outerCI.n, data=d) + theme_bw()

pdf(paste(fig.loc,title.f,".pdf", sep=""), width=w, height=h)
results$fig.results
dev.off()
m12 <- results[4]

stargazerAD(c(m11, m12), title="Effect on Approval, Hypothetical", filename=paste(fig.loc,title.f,"-table",".tex", sep=""))

smalld <- d[,c("asc", "asc0", "his", "pro", "prot", "com", "mob", "eli.f", "eli.c", 
               "authoritarian", "ally", "capabilities", "salience",  
               "pre.questions", "asc.or", "partner", 
               "na1.v", "na2.v", "na3.v", "na2.v.dn", "na3.v.dn", "gender.o", "educ.o", "age.o",
               "age.m", "gender.m", "educ.m")]

#, "start.time.n", "start.time.n2", "start.time.n3", "start.time.swd"
stargazer(smalld, title="Summary Statistics", out=paste(fig.loc, "summary-statistics-hypothetical",".tex", sep=""))




## Summary statistics for variables ##
asc ~ pro + prot + com + mob + eli.f + eli.c + 
  authoritarian + ally + capabilities + salience + his +
  pre.questions + asc.or +  partner + asc.or + asc0.v2  +
  na1.v2 + na2.v2 + na3.v2 + na2.v.dn + na3.v.dn + gender + educ + age + age.m + gender.m + educ.m + 
  start.time.n + start.time.n2 + start.time.n3 + start.time.swd 

pro.results <- rbind(pro.results, coef(summary(results$model))["pro",])
row.names(pro.results) <- c("approval", "approval, cov")

statement.results <- rbind(statement.results, coef(summary(results$model))[c("com", "mob", "his", "eli.f","eli.c"),])
l1 <- length(statement.results[,1])
row.names(statement.results) <- paste(row.names(statement.results), rep(", ", l1), c(rep("approval",l1/2), rep("approval, cov", l1/2)), sep="")




## Replicate analysis without inattentive respondents
## Added by RP on 9/18/2017
r <- !is.na(d$asc) & rb

type <- "Inattentive-Hyp"
outcome <- "Approval"
cont <- ""
title.w <- paste(type, ", Effect on ", outcome, cont, sep="")
title.f <- paste(type, "-", outcome, cont, sep="")

#coefficients.list<- c("his", "pro", "prot", "com", "mob", "eli.f", "eli.c")
coefficients.list<- c("his", "pro", "com", "mob", "eli.f", "eli.c")
model <- asc ~ pro + prot + com + mob + eli.f + eli.c + 
  authoritarian + ally + capabilities + salience + his +
  pre.questions + asc.or 
results <- results.fun(data=d[d$att==1,], r=r, model=model, title.w=title.w, fig.loc=fig.loc, coefficients=coefficients.list, innerCI.n=innerCI.n, outerCI.n=outerCI.n, w=w, h=h)
results[1]
results[2]

m11_inattentive <- results[4]

r <- !is.na(d$asc) & rb

type <- "Inattentive-Hyp"
outcome <- "Approval"
cont <- "Covariates"
title.w <- paste(type, ", Effect on ", outcome, ", Controlling ", cont, sep="")
title.f <- paste(type, "-", outcome, cont, sep="")

#coefficients.list<- c("his", "pro", "prot", "com", "mob", "eli.f", "eli.c")
coefficients.list<- c("his", "pro", "com", "mob", "eli.f", "eli.c")
model <- asc ~ pro + prot + com + mob + eli.f + eli.c + 
  authoritarian + ally + capabilities + salience + his +
  pre.questions + asc.or +  partner + asc.or + asc0.v2  +
  na1.v2 + na2.v2 + na3.v2 + na2.v.dn + na3.v.dn + gender + educ + age + age.m + gender.m + educ.m + 
  start.time.n + start.time.n2 + start.time.n3 + start.time.swd 
results <- results.fun(data=d[d$att==1,], r=r, model=model, title.w=title.w, fig.loc=fig.loc, coefficients=coefficients.list, innerCI.n=innerCI.n, outerCI.n=outerCI.n, w=w, h=h)
results[1]
results[2] 
results$fig.results

m12_inattentive <- results[4]

stargazerAD(c(m11_inattentive, m12_inattentive), title="Effect on Approval, Hypothetical (Inattentive respondents omitted)", filename=paste(fig.loc,title.f,"-table",".tex", sep=""))




##analysis of inattention distribution
r <- !is.na(d$asc) & rb

type <- "Random-Inattentive-Hyp"
outcome <- "Attentiveness"
cont <- ""
title.w <- paste(type, ", Effect on ", outcome, cont, sep="")
title.f <- paste(type, "-", outcome, cont, sep="")

#coefficients.list<- c("his", "pro", "prot", "com", "mob", "eli.f", "eli.c")
coefficients.list<- c("his", "pro", "com", "mob", "eli.f", "eli.c")
model <- att ~ pro + prot + com + mob + eli.f + eli.c + 
  authoritarian + ally + capabilities + salience + his +
  pre.questions + asc.or 
results <- results.fun(data=d, r=r, model=model, title.w=title.w, fig.loc=fig.loc, coefficients=coefficients.list, innerCI.n=innerCI.n, outerCI.n=outerCI.n, w=w, h=h)

m31_inattentive <- results[4]

r <- !is.na(d$asc) & rb

type <- "Random-Inattentive-Hyp"
outcome <- "Attentiveness"
cont <- "Covariates"
title.w <- paste(type, ", Effect on ", outcome, ", Controlling ", cont, sep="")
title.f <- paste(type, "-", outcome, cont, sep="")

#coefficients.list<- c("his", "pro", "prot", "com", "mob", "eli.f", "eli.c")
coefficients.list<- c("his", "pro", "com", "mob", "eli.f", "eli.c")
model <- att ~ pro + prot + com + mob + eli.f + eli.c + 
  authoritarian + ally + capabilities + salience + his +
  pre.questions + asc.or +  partner + asc.or + asc0.v2  +
  na1.v2 + na2.v2 + na3.v2 + na2.v.dn + na3.v.dn + gender + educ + age +
  start.time.n + start.time.n2 + start.time.n3 + start.time.swd 
results <- results.fun(data=d, r=r, model=model, title.w=title.w, fig.loc=fig.loc, coefficients=coefficients.list, innerCI.n=innerCI.n, outerCI.n=outerCI.n, w=w, h=h)

results[1]
results[2] 
results$fig.results


m32_inattentive <- results[4]

stargazerAD(c(m31_inattentive), title="Effect on Attentiveness, Hypothetical", filename=paste(fig.loc,title.f,"-table",".tex", sep=""))




#Checking if provocation or statement of commitment interacts with eli.f or eli.c.

model <- asc ~ pro + prot + com + mob + eli.f + eli.c + com*eli.f + com*eli.c + pro*eli.f + pro*eli.c +
  authoritarian + ally + capabilities + salience + his +
  pre.questions + asc.or +  partner + asc.or + asc0.v2  +
  na1.v2 + na2.v2 + na3.v2 + na2.v.dn + na3.v.dn + gender + educ + age + age.m + gender.m + educ.m + 
  start.time.n + start.time.n2 + start.time.n3 + start.time.swd 
coefplot(lm(model, data=d[r,]), title=title.w, 
         intercept=FALSE, lwdOuter=0.5, innerCI=innerCI.n, outerCI= outerCI.n, data=d) + theme_bw()

#Results: Same as above.


r <- !is.na(d$asc) & rb
title.w <- "Hyp, Effect on Resolve"
type <- "Hyp"
outcome <- "Resolve"
cont <- ""
title.w <- paste(type, ", Effect on ", outcome, cont, sep="")
title.f <- paste(type, "-", outcome, cont, sep="")

coefficients.list<- c("his", "pro", "com", "mob", "eli.f", "eli.c")
model <- rac ~ pro + prot + com + mob + eli.f + eli.c + 
  authoritarian + ally + capabilities + salience + his +
  pre.questions + rac.or
results <- results.fun(data=d, r=r, model=model, title.w=title.w, fig.loc=fig.loc, coefficients=coefficients.list, innerCI.n=innerCI.n, outerCI.n=outerCI.n, w=w, h=h)
results[1]
results[2] 
results$fig.results
coefplot(lm(model, data=d[r,]), title=title.w, 
         intercept=FALSE, lwdOuter=0.5, innerCI=innerCI.n, outerCI= outerCI.n, data=d) + theme_bw()

pdf(paste(fig.loc,title.f,".pdf", sep=""), width=w, height=h)
results$fig.results
dev.off()


#Interactions with elite cues?
model <- rac ~ pro + prot + com + mob + eli.f + eli.c + com*eli.f + com*eli.c + pro*eli.f + pro*eli.c +
  authoritarian + ally + capabilities + salience + his +
  pre.questions + rac.or
coefplot(lm(model, data=d[r,]), title=title.w, 
         intercept=FALSE, lwdOuter=0.5, innerCI=innerCI.n, outerCI= outerCI.n, data=d) + theme_bw()

#Provocation (maybe) positively interacts with elite cues for resolve. So provocation now has no main effect, but especially operates after elite cues. Why?  I don't know. Maybe the elite cues work when honor is not at stake. But when honor is at stake, they seem like wimply excuses. The positive interaction is stronger for elite cue costs. And similarly, a provocation can be more easily swallowed if you don't make excuses, but if you do it is perceived as worse. ?

#reduce resolve for those motivated by  


pro.results <- rbind(pro.results, coef(summary(results$model))["pro",])
row.names(pro.results)[length(pro.results[,1])] <- "resolve"

toprint <- paste("p_{+}=",round(pro.results[3,4]/2,3),  sep="")
write(toprint, file=paste(fig.loc,type,outcome,"p",".tex", sep=""))

statement.results <- rbind(statement.results, coef(summary(results$model))[c("com","mob", "his", "eli.f","eli.c"),])
ln <- length(statement.results[,1])-l1
r <- (l1+1):(l1+ln)
row.names(statement.results)[r] <- paste(row.names(statement.results)[r], rep(", ", ln), c(rep("resolve",ln)), sep="")



r <- !is.na(d$asc) & rb
title.w <- "Hyp, Effect on Resolve, Controlling Covariates"
type <- "Hyp"
outcome <- "Resolve"
cont <- "Covariates"
title.w <- paste(type, ", Effect on ", outcome, ", Controlling ", cont, sep="")
title.f <- paste(type, "-", outcome, cont, sep="")

#coefficients.list<- c("his", "pro", "prot", "com", "mob", "eli.f", "eli.c")
coefficients.list<- c("his", "pro", "com", "eli.f", "eli.c")
model <- rac ~ pro + prot + com + mob + eli.f + eli.c + 
  authoritarian + ally + capabilities + salience + his +
  pre.questions + rac.or +  partner + asc.or + asc0.v2  +
  na1.v2 + na2.v2 + na3.v2 + na2.v.dn + na3.v.dn + gender + educ + age + age.m + gender.m + educ.m + 
  start.time.n + start.time.n2 + start.time.n3 + start.time.swd 
results <- results.fun(data=d, r=r, model=model, title.w=title.w, fig.loc=fig.loc, coefficients=coefficients.list, innerCI.n=innerCI.n, outerCI.n=outerCI.n, w=w, h=h)
results[1]
results[2] 
results$fig.results
coefplot(lm(model, data=d[r,]), title=title.w, 
         intercept=FALSE, lwdOuter=0.5, innerCI=innerCI.n, outerCI= outerCI.n, data=d) + theme_bw()

pdf(paste(fig.loc,title.f,".pdf", sep=""), width=w, height=h)
results$fig.results
dev.off()


r <- !is.na(d$asc) & rb
#Interactions with elite cues?
model <- rac ~ pro + prot + com + mob + eli.f + eli.c + com*eli.f + com*eli.c + pro*eli.f + pro*eli.c +
  authoritarian + ally + capabilities + salience + his +
  pre.questions + rac.or +  partner + asc.or + asc0.v2  +
  na1.v2 + na2.v2 + na3.v2 + na2.v.dn + na3.v.dn + gender + educ + age + age.m + gender.m + educ.m + 
  start.time.n + start.time.n2 + start.time.n3 + start.time.swd 
coefplot(lm(model, data=d[r,]), title=title.w, 
         intercept=FALSE, lwdOuter=0.5, innerCI=innerCI.n, outerCI= outerCI.n, data=d) + theme_bw()
#Similar results. 



pro.results <- rbind(pro.results, coef(summary(results$model))["pro",])
row.names(pro.results)[length(pro.results[,1])] <- "resolve, cov"

statement.results <- rbind(statement.results, coef(summary(results$model))[c("com", "mob", "his", "eli.f","eli.c"),])
l1 <- ln+l1
ln <- length(statement.results[,1])-l1
r <- (l1+1):(l1+ln)
row.names(statement.results)[r] <- paste(row.names(statement.results)[r], rep(", ", ln), c(rep("resolve, cov",ln)), sep="")



r <- !is.na(d$asc) & rb
title.w <- "Hyp, Effect on Resolve2"
type <- "Hyp"
outcome <- "Resolve2"
cont <- ""
title.w <- paste(type, ", Effect on ", outcome, cont, sep="")
title.f <- paste(type, "-", outcome, cont, sep="")

#coefficients.list<- c("his", "pro", "prot", "com", "mob", "eli.f", "eli.c")
coefficients.list<- c("his", "pro", "com", "eli.f", "eli.c")
model <- ra3c ~ pro + prot + com + mob + eli.f + eli.c + 
  authoritarian + ally + capabilities + salience + his +
  pre.questions 
results <- results.fun(data=d, r=r, model=model, title.w=title.w, fig.loc=fig.loc, coefficients=coefficients.list, innerCI.n=innerCI.n, outerCI.n=outerCI.n, w=w, h=h)
results[1]
results[2] 
results$fig.results
coefplot(lm(model, data=d[r,]), title=title.w, 
         intercept=FALSE, lwdOuter=0.5, innerCI=innerCI.n, outerCI= outerCI.n, data=d) + theme_bw()

pdf(paste(fig.loc,title.f,".pdf", sep=""), width=w, height=h)
results$fig.results
dev.off()

r <- !is.na(d$asc) & rb
title.w <- "Hyp, Effect on Resolve2, Controlling Covariates"
type <- "Hyp"
outcome <- "Resolve2"
cont <- "Covariates"
title.w <- paste(type, ", Effect on ", outcome, ", Controlling ", cont, sep="")
title.f <- paste(type, "-", outcome, cont, sep="")

#coefficients.list<- c("his", "pro", "prot", "com", "mob", "eli.f", "eli.c")
coefficients.list<- c("his", "pro", "com", "eli.f", "eli.c")
model <- ra3c ~ pro + prot + com + mob + eli.f + eli.c + 
  authoritarian + ally + capabilities + salience + his +
  pre.questions +  partner + asc.or + asc0.v2  +
  na1.v2 + na2.v2 + na3.v2 + na2.v.dn + na3.v.dn + gender + educ + age + age.m + gender.m + educ.m + 
  start.time.n + start.time.n2 + start.time.n3 + start.time.swd 
results <- results.fun(data=d, r=r, model=model, title.w=title.w, fig.loc=fig.loc, coefficients=coefficients.list, innerCI.n=innerCI.n, outerCI.n=outerCI.n, w=w, h=h)
results[1]
results[2]
results$fig.results
coefplot(lm(model, data=d[r,]), title=title.w, 
         intercept=FALSE, lwdOuter=0.5, innerCI=innerCI.n, outerCI= outerCI.n, data=d) + theme_bw()

pdf(paste(fig.loc,title.f,".pdf", sep=""), width=w, height=h)
results$fig.results
dev.off()


pvalues <- -1


pdf(paste(fig.loc,"Hyp-Approval(Pre)-Histogram",".pdf", sep=""), width=w, height=h*1.6)
par(mar=c(3.1,4.1,4,2.1))
hist(d$asc0, main="Approval (Pre) in Hypothetical", breaks=seq(0.5, 5.5, 1), freq=FALSE, labels=TRUE, col="blue", xaxt='n', xlab="", ylim=c(0,0.5))
axis(side=1, at=c(1, 3, 5), labels=c("Strongly \n disapprove", "Neither approve \n nor disapprove", "Strongly \n approve"))
dev.off()

pdf(paste(fig.loc,"Hyp-Approval(Post)-Histogram",".pdf", sep=""), width=w, height=h*1.6)
par(mar=c(3.1,4.1,4,2.1))
hist(d$asc, main="Approval (Post) in Hypothetical", breaks=seq(0.5, 5.5, 1), freq=FALSE, labels=TRUE, col="blue", xaxt='n', xlab="", ylim=c(0,0.5))
axis(side=1, at=c(1, 3, 5), labels=c("Strongly \n disapprove", "Neither approve \n nor disapprove", "Strongly \n approve"))
dev.off()


pro.results <- as.data.frame(pro.results)
pro.results$names <- as.factor(row.names(pro.results))
pro.results$names2 <- as.factor(c("No Covariates", "Covariates", "No Covariates", "Covariates"))
pro.results$dv <- c(0,0,1,1)
names(pro.results) <- c("est", "se", "t", "p", "names", "names2", "dv")
pro.results

statement.results <- as.data.frame(statement.results)
statement.results$names <- as.factor(row.names(statement.results))
statement.results$names2 <- as.factor(rep(c(rep("No Covariates",5),rep("Covariates",5)),2))
statement.results$names3 <- as.factor(rep(c("Explicit threat", "Mobilization", "Nationalist History", "Biding Time", "Cost of War"),4))
statement.results$dv <- c(rep(0,10), rep(1,10))
names(statement.results) <- c("est", "se", "t", "p", "names", "names2","names3", "dv")
statement.results$names4 <- paste(statement.results$names3, "\n ", statement.results$names2, sep="")

statement.results

innerCI.n <- 1.64 #set so exclusion of CI implies two sided significance of 0.1
outerCI.n <- 1.96 #set so exclusion of CI implies two sided significance of 0.05

r <- !is.na(d$asc) & rb
n <- length(d[r,1])


type <- "Hyp"
r2 <- statement.results$dv==0 & (statement.results$names3=="Explicit threat" | statement.results$names3=="Mobilization")   #& statement.results$names2=="No Covariates" 
dv.n <- "Approval"
title.w <- "Effect on"
covs <- ""
main.t <- "1StateA-"
ymin.n <- -1; ymax.n <- .4
source("./6-coef-figuresb.R", echo=TRUE)

type <- "Hyp"
r2 <- statement.results$dv==0 & (statement.results$names3!="Explicit threat" & statement.results$names3!="Mobilization")   #& statement.results$names2=="No Covariates" 
dv.n <- "Approval"
title.w <- "Effect on"
covs <- ""
main.t <- "1StateB-"
ymin.n <- -1; ymax.n <- .4
source("./6-coef-figuresb.R", echo=TRUE)


type <- "Hyp"
r2 <- statement.results$dv==0 & statement.results$names2=="No Covariates"
dv.n <- "Approval"
title.w <- "Effect on"
covs <- ""
main.t <- "1State-"
ymin.n <- -1; ymax.n <- .4
source("./6-coef-figuresb.R", echo=TRUE)

type <- "Hyp"
r2 <- statement.results$dv==1 & statement.results$names2=="No Covariates"
dv.n <- "Resolve"
title.w <- "Effect on"
covs <- ""
main.t <- "1State-"
ymin.n <- -1; ymax.n <- .4
source("./6-coef-figuresb.R", echo=TRUE)

type <- "Hyp"
r2 <- statement.results$dv==0 & statement.results$names2=="Covariates"
dv.n <- "Approval"
title.w <- "Effect on"
covs <- "Covs"
main.t <- "1State-"
ymin.n <- -1; ymax.n <- .4
source("./6-coef-figuresb.R", echo=TRUE)

type <- "Hyp"
r2 <- statement.results$dv==1 & statement.results$names2=="Covariates"
dv.n <- "Resolve"
title.w <- "Effect on"
covs <- "Covs"
main.t <- "1State-"
ymin.n <- -1; ymax.n <- .4
source("./6-coef-figuresb.R", echo=TRUE)


type <- "Hyp"
r2 <- pro.results$dv==0
dv.n <- "Approval"
title.w <- "Effect of Provocation on"
main.t <- "1Prov-"
ymin.n <- -1; ymax.n <- .4
source("./6-coef-figures.R", echo=TRUE)

r2 <- pro.results$dv==1
dv.n <- "Resolve"
title.w <- "Effect of Provocation on"
ymin.n <- -1.5; ymax.n <- 3.2
source("./6-coef-figures.R", echo=TRUE)



#------------------------
#which country?
#----------------------


#which country? 
tab <- table(d$a.pcono)
tab2 <- prop.table(tab)
library(xtable)
xtable(tab2, digits=3)



#-----------------------------------------------------------------------
#National Honour 
#-------------------------------------------------------------------------


rm(list = ls())

setwd("~/Dropbox/Dafoe-Weiss-Empirics/18-04-12-ISQ-Raluca/Code/Replication File ISQ/")

# Customized Stargazer Function, without asterisk inflation
stargazerAD <- function(table, title="Default", filename="Figures/default.tex", dep.var.labels=NULL) {
  a <- stargazer(table, title=title,  out=filename, 
                 align=TRUE, star.cutoffs = c(0.1, 0.05, 0.01, 0.001), star.char = c("\\dagger", "*", "**", "***"), 
                 notes="$^{\\dagger} p<0.1$; $^{*} p<0.05$; $^{**} p<0.01$; $^{***} p<0.001$", notes.append=FALSE, single.row=TRUE)
  return(a)
}


load("./MASTER CODE/hypanalysis.Rdata")

#when prequestions asked?

#qplot(d$start.time,d$pre.questions, geom='smooth', span =0.9)  

ggplot(d,aes(x=start.time.n,y=pre.questions))+stat_smooth(se=F) + 
  geom_vline(xintercept = d$start.time.n[1850])  


df <- d[d$start.time.n<d$start.time.n[1850],]

m <-lm(asc~pre.questions, data=df)

stargazerAD(m, title="Prequestions Effect on Approval (Hypothetical Design)", file="../../Figures/hyp_preq.tex")



#-----------------------------------------------------------------------
#Raw Results
#-------------------------------------------------------------------------

rm(list = ls())
library(ggplot2)
library(coefplot)
library(stargazer)
library(scales)
library(tidyverse)
setwd("~/Dropbox/Dafoe-Weiss-Empirics/18-04-12-ISQ-Raluca/Code/Replication File ISQ/")

load("./MASTER CODE/hypanalysis.Rdata")

myvars <- names(d) %in% c("his", "pro", "prot", "com", "mob", "eli.f", "eli.c", "asc") 
small_data <- d[myvars]
small_data$ID <- seq.int(nrow(small_data))

#Constructing Appropriate Comparisons 
#Treatments are assigned independently, except only one of biding time or cost of war can occur 

# the control is those without biding time and those without cost of war (only one can happen)
small_data$new_eli.f[small_data$eli.f==1] <- 1
small_data$new_eli.f[small_data$eli.f==0 & small_data$eli.c==0] <- 0

#similarly for cost of war

small_data$new_eli.c[small_data$eli.c==1] <- 1
small_data$new_eli.c[small_data$eli.c==0 & small_data$eli.f==0] <- 0


#making a dataset for each treatment 

#history
histvars <- names(small_data) %in% c("his", "asc")
hist_data <- small_data[histvars]

hist_data$his_type[small_data$his==1] <- "Treatment"
hist_data$his_type[small_data$his==0] <- "Control"

sum(is.na(hist_data$his_type)) ## no missing for history treatment
sum(is.na(hist_data$asc)) ## some missing for asc. we'll take them out for the ggplot

hist_data <- na.omit(hist_data) ##taking out missing values 

hist_data2 <- hist_data %>% 
  group_by(his,asc) %>% 
  summarise(count=n()) %>% 
  mutate(perc=count/sum(count))

hist_data2$overall = "History"

hist_gg <-ggplot(hist_data2, aes(x=as.factor(asc), y=perc, group=as.factor(his), fill=as.factor(his)))
hist_gg2 <- hist_gg +
  geom_bar(stat= "identity", position="dodge") + 
  scale_y_continuous(labels=percent) + 
  facet_wrap(~overall)

hist_gg3 <- hist_gg2 +  
  ylab("Percentage") + 
  xlab (" ")  + 
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1))   + 
  scale_x_discrete(breaks=c(1,2,3,4,5),labels= c("S Disaprove", "Disaprove", "Neither", "Approve", "S Approve")) + 
  scale_fill_manual(labels=c("Control", "Treatment"), values=c("goldenrod1", "grey40"), name = "Treatment Status") + 
  coord_cartesian(ylim = c(0,.5))

#Commitment

comvars <- names(small_data) %in% c("com", "asc")
com_data <- small_data[comvars]

com_data$com_type[small_data$com==1] <- "Treatment"
com_data$com_type[small_data$com==0] <- "Control"

sum(is.na(com_data$com_type)) ## no missing 
sum(is.na(com_data$asc)) ## some missing for asc. we'll take them out for the ggplot

com_data <- na.omit(com_data) ##taking out missing values 

com_data2 <- com_data %>% 
  group_by(com,asc) %>% 
  summarise(count=n()) %>% 
  mutate(perc=count/sum(count))

com_data2$overall = "Commitment"

com_gg <-ggplot(com_data2, aes(x=as.factor(asc), y=perc, group=as.factor(com), fill=as.factor(com)))
com_gg2 <- com_gg +
  geom_bar(stat="identity", position="dodge") + 
  scale_y_continuous(labels=percent) + 
  facet_wrap(~overall)

com_gg3 <- com_gg2 +  
  ylab(" ") + 
  xlab (" ")  + 
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1))   + 
  scale_x_discrete(breaks=c(1,2,3,4,5),labels= c("S Disaprove", "Disaprove", "Neither", "Approve", "S Approve")) + 
  scale_fill_manual(labels=c("Control", "Treatment"), values=c("goldenrod1", "grey40"), name = "Treatment Status") + 
  coord_cartesian(ylim = c(0,.5))


#Mobilization

mobvars <- names(small_data) %in% c("mob", "asc")
mob_data <- small_data[mobvars]

mob_data$mob_type[small_data$mob==1] <- "Treatment"
mob_data$mob_type[small_data$mob==0] <- "Control"

sum(is.na(mob_data$mob_type)) ## no missing 
sum(is.na(mob_data$asc)) ## some missing for asc. we'll take them out for the ggplot

mob_data <- na.omit(mob_data) ##taking out missing values 
mob_data2 <- mob_data %>% 
  group_by(mob,asc) %>% 
  summarise(count=n()) %>% 
  mutate(perc=count/sum(count))


mob_data2$overall = "Mobilization"

mob_gg <-ggplot(mob_data2, aes(x=as.factor(asc), y=perc, group=as.factor(mob), fill=as.factor(mob)))
mob_gg2 <- mob_gg +
  geom_bar(stat="identity", position="dodge") + 
  scale_y_continuous(labels=percent) + 
  facet_wrap(~overall)

mob_gg3 <- mob_gg2 +  
  ylab(" ") + 
  xlab (" ")  + 
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1))   + 
  scale_x_discrete(breaks=c(1,2,3,4,5),labels= c("S Disaprove", "Disaprove", "Neither", "Approve", "S Approve")) + 
  scale_fill_manual(labels=c("Control", "Treatment"), values=c("goldenrod1", "grey40"), name = "Treatment Status") + 
  coord_cartesian(ylim = c(0,.5))


#Biding Time (eli.f)

elifvars <- names(small_data) %in% c("new_eli.f", "asc")
elif_data <- small_data[elifvars]

elif_data$elif_type[small_data$new_eli.f==1] <- "Treatment"
elif_data$elif_type[small_data$new_eli.f==0] <- "Control"

sum(is.na(elif_data$elif_type)) ## no missing 
sum(is.na(elif_data$asc)) ## some missing for asc. we'll take them out for the ggplot

elif_data <- na.omit(elif_data) ##taking out missing values 

elif_data2 <- elif_data %>% 
  group_by(new_eli.f,asc) %>% 
  summarise(count=n()) %>% 
  mutate(perc=count/sum(count))

elif_data2$overall = "Biding Time"

elif_gg <-ggplot(elif_data2, aes(x=as.factor(asc), y=perc, group=as.factor(new_eli.f), fill=as.factor(new_eli.f)))
elif_gg2 <- elif_gg +
  geom_bar(stat="identity", position="dodge") + 
  scale_y_continuous(labels=percent) + 
  facet_wrap(~overall)

elif_gg3 <- elif_gg2 +  
  ylab("Percentage") + 
  xlab (" ")  + 
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1))   + 
  scale_x_discrete(breaks=c(1,2,3,4,5),labels= c("S Disaprove", "Disaprove", "Neither", "Approve", "S Approve")) + 
  scale_fill_manual(labels=c("Control", "Treatment"), values=c("goldenrod1", "grey40"), name = "Treatment Status") + 
  coord_cartesian(ylim = c(0,.5))

#Cost of war (eli.c)

elicvars <- names(small_data) %in% c("new_eli.c", "asc")
elic_data <- small_data[elicvars]

elic_data$elic_type[small_data$new_eli.c==1] <- "Treatment"
elic_data$elic_type[small_data$new_eli.c==0] <- "Control"

sum(is.na(elic_data$elic_type)) ## no missing 
sum(is.na(elic_data$asc)) ## some missing for asc. we'll take them out for the ggplot

elic_data <- na.omit(elic_data) ##taking out missing values 

elic_data2 <- elic_data %>% 
  group_by(new_eli.c,asc) %>% 
  summarise(count=n()) %>% 
  mutate(perc=count/sum(count))

elic_data2$overall = "Cost of War"

elic_gg <-ggplot(elic_data2, aes(x=as.factor(asc), y=perc, group=as.factor(new_eli.c), fill=as.factor(new_eli.c)))
elic_gg2 <- elic_gg +
  geom_bar(stat="identity", position="dodge") + 
  scale_y_continuous(labels=percent) + 
  facet_wrap(~overall)

elic_gg3 <- elic_gg2 +  
  ylab(" ") + 
  xlab ("Government Approval")  + 
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1))   + 
  scale_x_discrete(breaks=c(1,2,3,4,5),labels= c("S Disaprove", "Disaprove", "Neither", "Approve", "S Approve")) + 
  scale_fill_manual(labels=c("Control", "Treatment"), values=c("goldenrod1", "grey40"), name = "Treatment Status") + 
  coord_cartesian(ylim = c(0,.5))


#putting them together

library(ggpubr)

ggarrange(hist_gg3, com_gg3, mob_gg3, elif_gg3, elic_gg3, 
          ncol = 3, nrow = 2,common.legend = TRUE, legend = "right", align = "v")



ggsave(file="../../Figures/hyp_vertical_only5.eps")


#--------------------------------------------------------------------------------------- 
#Coefficients Plot from Regression

#Regression Plot

treat1 <- lm(asc~his, data=hist_data)
treat4 <- lm(asc~com, data=com_data)
treat5 <- lm(asc~mob, data=mob_data)
treat6 <- lm(asc~new_eli.c, data=elic_data)
treat7 <- lm(asc~new_eli.f, data=elif_data)


mlplot <- multiplot(treat1,  treat4, treat5, treat6, treat7,
                    predictors = c("his",  "com", "mob", "new_eli.c","new_eli.f"), title="Effect on Approval, Hyp", xlab="Change in Approval")

mlplot + theme_bw() + 
  scale_y_discrete(label=c("History", "Commitment", 
                           "Mobilization", "Cost of War", "Biding Time")) + 
  theme(legend.position="none") + 
  scale_color_manual(values=c("black", "black", "black", "black", "black", "black", "black"))
ggsave(file="../../Figures/hyp_raw_coefplot_only5.eps")


#Error Bar Plot

#History

gg1 <- ggplot(hist_data, aes(as.factor(his_type), asc, fill=as.factor(his_type)))+ 
  stat_summary(fun.y=mean, geom="bar", width=1, position= "dodge") +
  stat_summary(fun.data =mean_se, geom="errorbar",  width=.2, size=1, color="black",
               position=position_dodge(width = .7)) +
  ylab(" ") +   xlab("History") + 
  scale_fill_manual(labels=c("Control", "Treatment"), values=c("goldenrod1", "grey40"), name="Treatment Status") + 
  coord_cartesian(ylim=c(2.7,3.2)) + 
  theme(panel.background = element_blank(), 
        axis.line.y = element_line(colour = "Black", linetype = "solid"),
        axis.ticks.length=unit(0.3,"cm"),
        axis.ticks.x = element_blank(),
        axis.text.x = element_blank(), 
        axis.line.x = element_blank()) +
  scale_y_continuous(limits = c(2.7,3.2),  expand = c(0,0), oob=rescale_none) 


#Commitment

gg4 <- ggplot(com_data, aes(as.factor(com_type), asc, fill=as.factor(com_type)))+ 
  stat_summary(fun.y=mean, geom="bar", width=1, position= "dodge") +
  stat_summary(fun.data =mean_se, geom="errorbar",  width=.2, size=1, color="black",
               position=position_dodge(width = .7)) +
  ylab(" ") +   xlab("Commit") + 
  scale_fill_manual(labels=c("Control", "Treatment"), values=c("goldenrod1", "grey40"), name="Treatment Status") + 
  scale_y_continuous(limits = c(2.7,3.2),  expand = c(0,0), oob=rescale_none)  + 
  coord_cartesian(ylim=c(2.7,3.2)) + 
  theme(panel.background = element_blank(), 
        axis.line.y = element_blank(),
        axis.ticks.length=unit(0.3,"cm"),
        axis.ticks.x = element_blank(),
        axis.text.x = element_blank(), 
        axis.ticks.y =element_blank(), 
        axis.text.y = element_blank(),
        axis.line.x = element_blank())

#Mobilization

gg5 <- ggplot(mob_data, aes(as.factor(mob_type), asc, fill=as.factor(mob_type)))+ 
  stat_summary(fun.y=mean, geom="bar", width=1, position= "dodge") +
  stat_summary(fun.data =mean_se, geom="errorbar",  width=.2, size=1, color="black",
               position=position_dodge(width = .7)) +
  ylab(" ") +   xlab("Mobilize") + 
  scale_fill_manual(labels=c("Control", "Treatment"), values=c("goldenrod1", "grey40"), name="Treatment Status") + 
  scale_y_continuous(limits = c(2.7,3.2),  expand = c(0,0), oob=rescale_none)  + 
  coord_cartesian(ylim=c(2.7,3.2)) + 
  theme(panel.background = element_blank(), 
        axis.line.y = element_blank(),
        axis.ticks.length=unit(0.3,"cm"),
        axis.ticks.x = element_blank(),
        axis.text.x = element_blank(), 
        axis.ticks.y =element_blank(), 
        axis.text.y = element_blank(),
        axis.line.x = element_blank())


#Biding Time

gg6 <- ggplot(elif_data, aes(as.factor(elif_type), asc, fill=as.factor(elif_type)))+ 
  stat_summary(fun.y=mean, geom="bar", width=1, position= "dodge") +
  stat_summary(fun.data =mean_se, geom="errorbar",  width=.2, size=1, color="black",
               position=position_dodge(width = .7)) +
  ylab(" ") +   xlab("Biding") + 
  scale_fill_manual(labels=c("Control", "Treatment"), values=c("goldenrod1", "grey40"), name="Treatment Status") + 
  scale_y_continuous(limits = c(2.7,3.2),  expand = c(0,0), oob=rescale_none)  + 
  coord_cartesian(ylim=c(2.7,3.2)) + 
  theme(panel.background = element_blank(), 
        axis.line.y = element_blank(),
        axis.ticks.length=unit(0.3,"cm"),
        axis.ticks.x = element_blank(),
        axis.text.x = element_blank(), 
        axis.ticks.y =element_blank(), 
        axis.text.y = element_blank(),
        axis.line.x = element_blank())

#Cost of War 

gg7 <- ggplot(elic_data, aes(as.factor(elic_type), asc, fill=as.factor(elic_type)))+ 
  stat_summary(fun.y=mean, geom="bar", width=1, position= "dodge") +
  stat_summary(fun.data =mean_se, geom="errorbar",  width=.2, size=1, color="black",
               position=position_dodge(width = .7)) +
  ylab(" ") +   xlab("Cost of War") + 
  scale_fill_manual(labels=c("Control", "Treatment"), values=c("goldenrod1", "grey40"), name="Treatment Status") + 
  scale_y_continuous(limits = c(2.7,3.2),  expand = c(0,0), oob=rescale_none)  + 
  coord_cartesian(ylim=c(2.7,3.2)) + 
  theme(panel.background = element_blank(), 
        axis.line.y = element_blank(),
        axis.ticks.length=unit(0.3,"cm"),
        axis.ticks.x = element_blank(),
        axis.text.x = element_blank(), 
        axis.ticks.y =element_blank(), 
        axis.text.y = element_blank(),
        axis.line.x = element_blank())



ggarrange(gg1, gg4, gg5, gg6, gg7,
          ncol = 5, nrow = 1,common.legend = TRUE, legend = "right")

ggsave(file="../../Figures/hyp_subst_only5.eps")


#End.
